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Abstract. The authors previously considered a method solving optimization 
problems by using a system of interconnected network of two component Bose-Einstein 
condensates (Byrnes, Yan, Yamamoto New J. Phys. 13, 113025 (2011)). The use of 
bosonic particles was found to give a reduced time proportional to the number of bosons 
N for solving Ising model Hamiltonians by taking advantage of enhanced bosonic 
cooling rates. In this paper we consider the same system in terms of neural networks. 
We find that up to the accelerated cooling of the bosons the previously proposed 
system is equivalent to a stochastic continuous Hopfield network. This makes it clear 
that the BEG network is a physical realization of a simulated annealing algorithm, with 
an additional speedup due to bosonic enhancement. We discuss the BEG network in 
terms of typical neural network tasks such as learning and pattern recognition and find 
that the latter process may be accelerated by a factor of N. 
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1. Introduction 

With tlie exception of a small class of problems that are solvable analytically, most 
quantum many problems can only be examined using numerical means, for which exact 
simulations scale exponentially with the problem size. Approximate methods such 
as quantum Monte Carlo and Density Matrix Renomalization Group (DMRG) give 
accurate results for certain cases but no general algorithm exists that can be applied 
to an arbitrary system. In the field of quantum information technology, quantum 
simulation has gathered a large amount of attention as an alternative means to study 
such quantum many body problems. A quantum simulator is a device where a quantum 
many body problem of interest is artificially created in the lab, such that its properties 
can be controlled and measured directly [1]. Via directly using quantum mechanics 
in the simulation, there is no exponential overhead in keeping track of the number of 
states in the Hilbert space of the problem. This is in the spirit of Feynman's original 
motivations for quantum computing [2] , where quantum mechanics, rather than classical 
mechanics, is used to simulate quantum many body problems. 

Given this general approach to quantum many problems, the question of whether 
a quantum simulation approach can be applied to a general Ising model becomes an 
important question. The Ising model problem consists of finding the lowest energy state 
of the Hamiltonian 



where Jij is a real symmetric matrix that specifies the connections between the sites 
and cTj = ±1 is a spin variable. The task is then to find the minimal energy spin 
configuration {cTj} of the Hamiltonian ([T]). The problem of finding the solution of the 
Hamiltonian ([1]) is in fact known to be NP-complete, since it can be trivially be mapped 
onto the MAX-CUT problem [3]. Furthermore, it can in principle encode an arbitrary 
NP-complete problem by a polynomial time mapping procedure thus the potential 
application of a reliable method of solving ([1]) is extremely broad. Although ([1]) is 
itself a classical Hamiltonian since it does not contain any non-commuting operators, 
as with quantum annealing where the Hamiltonian is supplemented with an additional 
transverse field, quantum "tricks" may be used to speed up the solution of the ground 
state beyond classical methods. 

In a previous work we investigated a computational device which finds the solution 
of an Ising model by a set of interconnected Bose-Einstein condensates [3, E]. In the 
approach of Ref. |5], each spin was replaced by a system of bosons which can take one 
of two states. By implementing an analogous Hamiltonian to (II]) and cooling the system 
down into into ground state, it was shown that the solution of the original Ising model 
problem could be found. There is a speedup compared to simply implementing ^ using 
single spins, because of the presence of bosonic final state stimulation within each bosonic 
spin. This resulted in finding the solution of ^ with a speedup of approximately A^. 
The attractive feature of the proposal in Ref. [5] is that the computation is done simply 




(1) 



i,j=l i=l 
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by implementing a static Hamiltonian, so that no complicated gate sequence needs to 
be employed in order to find the ground state solution. Effectively, the dissipative 
process itself performs the computation itself, and therefore can also be considered to 
be a reservoir engineering approach to quantum computation. Related proposals were 
offered in Refs. [TJ IH], where instead of BECs, photons were used. 

In this paper, we analyze the proposal in Ref. [5] from the point of view of neural 
networks, specifically the stochastic continuous Hopfield model. Recasting the proposal 
in this form allows for a clearer analysis of the properties of the device, where standard 
results can be carried over to the BEC case. It clarifies the origin of the ~ N speedup of 
the device, which was established via a numerical approach in Ref. [5]. We find that the 
~ speedup originates from each element of the Hopfield network being accelerated 
due to bosonic stimulation, and thermal fluctuations provide the stochastic aspect to 
the Hopfield network. We then consider some simple applications of the BEC network 
for neural networking tasks. 

2. Bose-Einstein condensate networks 

In order to make this paper self-contained, we first give a brief description of the proposal 
of Ref. [5] . In order to solve ^ we consider a system such as that shown in Figure [1] 
Each spin cxj in Hp is associated with a trapping site containing bosonic particles. 
The bosons can occupy one of two spin states, which we label by a = ±1. 



where is the annihilation operator for spin a on site i. On each site i we may define 
a total spin operator Si = a|^ai+ — a|_ai_ taking eigenvalues Si\k) = {—N + 2k)\k). 
The sites are controlled such that the system follows the bosonic Ising Hamiltonian 



where Jij is the same matrix as in Hp which specifies the computational problem. 
One such implementation of creating the interactions is via a measurement-feedback 
approach. In this approach, the total spin Si on each site is continuously measured, 
and an appropriate control signal is continuously fed back into the system by applying 
a local dc field on another site. Given a measurement result of {Sj(t)} across the spins, 
a local field 




(2) 



M M 

H = ^ JijSiSj + ^ XiSi 

i,j=l i=l 



(3) 



B,{t)=J2JijSjit) + X, 



(4) 



j 

is applied on site i. Since the Zeeman energy due to this field is 



H = ^ BiSi, 



(5) 



a simple substitution yields ([3]). 
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Figure 1. Each site of tlie Ising Hamiltonian is encoded as a trapping site, containing 
N bosons. The bosons can occupy one of two states tr = ±1, depicted as either red 
or blue. The interaction between the sites may be externally induced by measuring 
the total spin Si on each site i via the detectors. A local field on each site equal to 
Bi ~ JijSj + Ai is applied via the feedback circuit. 



Starting from a random spin configuration {S'j(t)}, the system is cooled assuming 
that the ambient temperature T is fixed. The procedure is essentially identical to a 
simulated annealing procedure, the sole difference being the use of the bosonic Ising 
model. By varying the tempeature during the cooling process such that it is time 
dependent T(t), a strategy similar to thermal annealing may be performed, in order to 
escape being trapped in local minima. In practice, instead of varying the temperature, 
varying the overall magnitude of ([3]) by adjusting the strength of the magnetic field (j4]) is 
equivalent. By taking advantange of the bosonic amplification of the cooling process, it 
was found in Ref. |5] that an approximate factor of N was found in the cooling process. 

The time evolution is modeled by a extension of the method presented by Glauber 
(ref. 14) to bosons. Given the M site Hamiltonian ([3]), the states are labeled 

M 1 

\k) = n ^— — (4+)'^(4-)"^-''l0), (6) 



=iy^A;,!(iV-A;,)! 



where the fcj range from to N , a\„ is the creation operator for a boson on site i in the 
state 0", and we have defined the vector k = {ki, ^2, • • • , /^m)- At t = 0, there is an equal 
probability of all states \k). The probability distribution then evolves according to 



dpk 
dt 



M 



J2 -w{k, di)pk+w{k+5i, -5i)pk+s,-w{k, -Si)pk+w{k-di, 5i)pk-s, 

i=l 

where di is a vector of length M with its ith element equal to one and zero otherwise. The 
w{k, di) is a weight factor representing the transition |fc) |fc + <5j), containing both the 
bosonic final state stimulation factor and a coefficient to ensure that the system evolves 
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to the correct thermal equihbrium distribution. We have restricted the transitions to 
first order transitions in ([7j) for simphcity. The final state stimulation factor can be 
calculated by assuming a Hamiltonian Ht = J2ia oLaj-o- and calculating the transition 
rate according to Fermi's golden rule (up to a constant) 

\{k + d,\Ht\k)\^ = {h + l){N-h), (8) 

At thermal equilibrium, the transition rates are equal between the states |fc) ^ \k-\-di), 
which ensures that ^ = 0. The final state stimulation factors cancel for such pairs and 
we have 

w(fc,(5,) ^ pk+s, . 
w{k + Si,-Si) Pk 

and similarly for Si — )■ — From the probability distribution at thermal equilibrium 
we can calculate 



Pk+Sj 

Pk 



exp 



-2/3 



Y, Jij{N -2k,) + Xi 



(10) 



where (3 = This gives the coefficients as 



w{k,di) =a(l + 7i)(A;. + l)(iV-A;,) 



:iii 



where 



7j = tanh —(5 



Y,J,^{N-2k,) + \ , (12) 

J^i J / 

a is a constant determining the overall transition time scale, a is set to 1 for all numerical 
calculations. 



3. Equivalence to the stochastic continuous Hopfield model 

We now show that the scheme detailed in the previous section is formally equivalent to 
the stochastic continuous Hopfield model. After defining the Hopfield model we show 
the equivalences between the two systems, and derive the evolution equations for the 
BEC network in the context of the Hopfield model. 



3.1. Definition of the Hopfield model 

The Hopfield model is an asynchronous model of an artificial neural network [9l [TO] , 
where the each unit of the network changes state at random times, such that no two units 
change simultaneously. The Hopfield model is in the class of recurrent neural networks, 
where the output of the network is fed back into itself, so that given a particular initial 
configuration, the system undergoes a sequence of changes in configuration of the units 
until steady state is achieved. 
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In the standard (discrete, deterministic) Hopfield model, each unit takes the value 
(jj = ±1. The units are updated according to the rule 



(Ti ^ = sgn 



(13) 



where sgn(x) = x/\x\ denotes the sign function, Wij is a symmetric matrix with zero 
diagonal elements, hi is the threshold of unit i. The units are updated to their new 
values cr[ at random. Whether a given state is stable with respect to the updates can 
be determined by the Lyapunov (energy) function 

^ = E W,,ai<jj + E h<yi. (14) 

ij i 

The update sequence proceeds until a local minima with respect to the Lyanupov 
function is reached. From a physics perspective, the transition rule (JT3ll can be 
viewed as a cooling process at zero temperature, where given an initial high energy 
spin configuration, the spins cool one by one randomly into a low energy state. It is 
thus equivalent to the discrete Ising model problem ([1]) up to a contant energy factor 
originating from the diagonal elements of Jij , which play no part in the dynamics of the 
problem. 

The model can be extended to one where continuous variables Xi = [—1, 1] are 
used in place of the discrete ones cTj. Such continuous Hopfield networks have similar 
properties to the discrete version in terms of the configuration of stable states [11]. The 
way the model is usually defined is in terms of an electric circuit model (see for e.g. Ref. 
[To]). On a single unit i, the time evolution of the circuit obeys 

* j 

where Ci is the capacitance, Ri is the resistance, Vi is the voltage, and Xi is the output 
of the circuit after operation of the nonlinear amplifier (or activation function). The 
activation function restricts the output to a limited region such that output is always 
in the range [—1, 1]. A typical choice is 

Xi = (pivi) = tanh(t>j). (16) 

The corresponding energy function for the dynamics is 

E = Y, W,,x,x, + ^ hx, + J2^.J <P~\x)dx. (17) 

ij i j * 

From the energy function it may then be shown that the system is guaranteed to converge 
to a local minima of the energy. Due to the energy structure of the continuous model 
being equivalent to the discrete model [in], a solution to the continuous model then 
gives a one-to-one correspondence to the discrete model. 
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3.2. Equations of motion on a single site 

In order to see the equivalence with the Hopfield model, let us first examine the dynamics 
on a single site, and set M = 1. In this case the probability distribution of the states 
evolve as 

^ = -w{k, l)pk + w{k + 1, -l)pk+i - w{k, -l)pk + w{k - 1, l)pk-i- 
dt 



where 

w{k, 1) = a(l + 7)(A; + 1){N - k) 
w{k, -1) = a(l - 7)A;(A^ - k + 1) 
7 = tanh (— /3A) , 

in this case since there is only one site. Multiplying the whole equation by k and 
summing over k given an equation for the mean value 

M = a(i + + l)(iV - k)) - ail - ^){k{N -k + 1)). (19) 
dt 

Making the approximation that (A;^) ~ (k)'^, and changing variables to S* = — + 2k 
the equation can be recast into the form 
1 d^ 

= -Njs^ -2s + -f{2 + N) (20) 

a dt 

where we have used the normalized variable s = {S)/N. An explicit solution for this 
may be found to be 

s{t) = A tanh {a-fANt + Kq) - ^ (21) 



where A = y 1 + ]v + jvV ^ constant that is order unity at low temperatures such 
that A^7 ^ 1. Kq is fixed by the initial conditions, for which we typically assume 
s(t = 0) = 0. We see that the spin approaches its steady state value with a timescale 
~ l/a|7|A^. At zero temperature where 7 = ±1 (depending upon the direction of the 
applied external field), the time scale is enhanced by a factor of A^, which is due to 
bosonic final state stimulation, as found in Ref. [5|. 

Let us compare this to the Hopfield network for a single site. In this case the 
equation of motion reads 

Changing variables to the output of the nonlinear amplifier, we obtain 



C 



dx ( dd) ^\ ^ 



dt \ dx 



(23) 



R 

This is the counterpart of the equation fl2Ul) for the Hopfield case. The solution of this 
x{t) = (f){bR + exp[-t/CR + K'q]). (24) 
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where K'q is set by the initial conditions. 

Writing the equations in this form makes the correspondence between the BEC 
system and the Hopfield network clear, which we summarize in Table [H The output x 
of the Hopfield network corresponds to the normalized spin variable s. The magnetic 
field B applied on each site for the BEC system then corresponds to the voltage on each 
Hopfield unit before the non-linear amplifier. The steady state values are determined 
by setting the right hand side of fl20|) and fl23l) to zero. The value of the steady state 
depends upon the ratio of the field (A or h) applied and the temperature ksT = or 
the conductance 1/R respectively. The overall time constant is controlled by 1/a or C in 
each case. In the Hopfield network, there is obviously no bosonic final state stimulation, 
so the speedup proportional to is absent. While the exact equation of motion for the 
two cases ( !20|) and ( l23l) differ in their details, it is clear that the qualitatively there is 
a similar structure and behavior to the dynamics. While in the Hopfield network, the 
overall time constant is determined by the capacitance of the circuit, the fundamental 
timescale of the BEC network is determined by the cooling rate of the bosons on each 
site. 

The analogue of the activation function 0(f) may be derived as follows. Considering 
the activation function as a rule that converts the internal magnetic field B to the spin 
variable s, we may derive the average spin at thermal equilibrium using the Boltzmann 
distribution taking into account of bosonic statistics. From the partition function 

pl^ = (1 - exp(-2/35)) exp(-2/35A;), (25) 

and therefore 



32-(2+^)) + (2 + Ar)(e2-(^+i) 



e 



2z\ 



= ^(-) (l-e2^)(l-e2-a^^)) ^''^ 

where z = BP here. The function above has a dependence that has similar behavior to 

^z) ^ tanh (^- '^^^'^^ ^ (27) 

which makes it clear that it plays a similar role to that of the activation function in 
the Hopfield network. 

3. 3. Equations of motion for interconnected BECs 

We may now generalize to the multi-site case. Multiplying by ki and summing over 
all k gives the equations 

^ = «((1 + 7.(fc))(fci + l)(iV - k,)) - a((l - 7.(fc))A;,(iV -k + 1)).(28) 

where we have written 7^ — )■ 7i(fc) to remind ourselves that this is not a constant in this 
case. Making the approximation that (fc™") ^ (A;)™", and making the change of variables 
to Si = {Si)/N, we obtain 
1 ds- 

--^ = -N^^is)s^ - 2s, + %is){2 + N). (29) 
a at 
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Quantity 



BEC network Hopfield model 



Steady-state control parameter 



Activation function 



Site field variable 
Ising matrix 
Local field 
Time constant 



Site variable 




Table 1. Equivalences between the network of BECs proposed in Ref. [5| and the 
Hopfield model. 



The sole difference to the single site case here is that the equilibrium values of Si are 
now dependent on the spins of all the other sites s. The dynamics on each site is the 
same as the single site case, and thus evolves in time as fl^T]) . considering the other spins 
to be approximately fixed. This basic structure is precisely the same dynamics that 
determine the equation of motion of the Hopfield network (|T5l) . 

Although it is not possible to solve the set of equations fl29|) analytically, we may 
see in this formulation why the whole system should have a speedup of ~ A^, as found 
in Ref. [5]. Considering a asynchronous update procedure (in fact this is exactly what 
was performed to simulate the dynamics of the system in the Monte Carlo procedure 
j6]), then all but one "active" site is fixed in spins. The active site then evolves in time 
according to the evolution of (|2T|) . This has a speedup of ~ A^|7| in the evolution of 
the spin, thus to make an incremental change 5s in the active spin takes a time reduced 
by A^|7| compared to the A^ = 1 case. The spin is then fixed, and then another site is 
chosen at random and this is updated. Since each step takes a reduced time of A^|7|, 
the whole evolution proceeds at an aceelerated rate. For sufficiently low temperatures, 
7 ^ ±1, and therefore the speedup is approximately ~ A^. 

3.4- Stochastic Hopfield network 

Up to this point, the equations of motion ( 129|) have been entirely deterministic. The role 
of the temperature was to merely shift the equilibrium values of the spins, as determined 
by (!26l) . and did not contribute to any stochastic evolution of the system. In the BEC 
network there are thermal fluctuations which cause the system to behave in a random 
manner. Therefore in order to fully capture the dynamics of the BEC network we 
must include the contribution of the random thermal process. Such stochastic versions 
of Hopfield networks, and their generalization to the Boltzmann machine (by having 
additional hidden units) are defined by modifying the update rule f|T3|) to include 
probabilistic effects. Considering the discrete Hopfield model first, the algorithm consists 
of selecting a particular spin cxj, and making an update according to 



if < 
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(Ji -> -(Ti if > with probability c'^'^/^bT 

o'i if AE > and otherwise. (30) 

AE is the energy difference between the state with the flipped spin and no flipped spin 

m- 

The evolution equations ([7]) for the BEC network can be converted into a set of 
stochastic update rules which give the same time depedence when an ensemble average 
is taken. The stochastic formulation also allows for a convenient method of numerically 
simulating the system, which was discussed in detail in Ref . |6] . We briefly describe the 
procedure as applied to the current formulation of the BEC network. The simulation 
is started from a random initial value of = {ki, k2, /ca/) in (jZ]), and we update the 
system by repeating the stochastic transition process following the kinetic Monte Carlo 
method [13]. In each update we calculate the transition weight w{k,di) in (JTj) for all 
the possible transitions. The transition is then made with a probability in proportion 
to the transition weight w{k, 6i). The time increment is then calculated according to 

At = -ln(r)M„, (31) 

where r G (0, 1] is a randomly generated number and 

i 

This procedure is repeated for many trajectories so that ensemble averages of quantities 
such as the average spin can be taken. 

This procedure produces exactly the same time dynamics as ([7]), hence for quantities 
such as the equilibration time this procedure must be followed. However, if only 
the behavior at equilibrium is required, the update procedure can be replaced by the 
Metropolis algorithm. The update procedure is then as follows. Start from a random 
initial value of fc = {ki, k2, A^Af). Then make an update according to 

ki^ ki±l if < 

ki-^ ki±l if AE > with probability e'^^/''^^ 

ki — )■ ki if AE > and otherwise. (33) 

where the energy difference in this case is 



AE = ±2 



(34) 



Y^J,^(N-2k,) + X, 

The exponential factor is precisely the same as that determining the weight factors in 
( ITOj) . The only difference in this case is that the bosonic stimulation factors are not 
present in this case, which are important only for determining the transition rates, and 
not the equilibrium values. The above Metropolis transition rule ffTOj) is identical to the 
stochastic Hopfield network, up to the difference that each site contains energy levels 
between ki = 0, . . . , N . It is therefore evidents in this context that the two systems are 
equivalent in their dynamics. 
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We may also derive the equivalent stochastic differential equations for the average 
spin by adding noise terms to fl29|) . First consider a single site, and start with the 
probability distribution ( |T8l) . Assuming that ^ 1, introduce the variables z = k/N, 
e = 1/N, and the density q{z,t) = pk/e of z at time t, the master equation is rewritten 

dq{z,t) 



dt 



-w{z, e)q{z, t)+w{z+e, -e)q{z+e, t)-w{z, -e)q{z, t)+w{z-e),t){35) 



Expanding w{z ± e, =Fe) and q{z ± e, t) up to second order in e, we obtain 



dq{z,t) 
dt 



d 

- [{w{z, e) - w{z, e))eq{z, t)] 



+ 



1 



{w{z, e) + w{z, -e))e^q{z, t)] + 0{e^) 
1 



2&2 

-^[A{z)q{z,t)] 



2dz^ 



[B,{z)q{z,t)]+0{e' 



(36) 



where 



A, {z) = {w{z, e) - w{z, -e))e 

B, iz) = iwiz,e)+w{z,-e))e^. (37) 

Using the diffusion approximation such that the transition rates are on the order of 
w{z,±e) ~ 1/dt, {w{z,e) ~w{z,e)) ~ e/dt, and e'^/dt ~ 0(1), and taking the limits of 
e — > 0, we obtain the Fokker-Planck equation 



dq(z,t) d , , 19^ , , 



(3S 



where A{z) = lime^o^e(^) and B{z) 
differential equation is given by 

dz 

It 



A{z) + jB{z)m 



\im.^^QB^{z). The corresponding stochastic 



(39) 



where ^{t) is Gaussian white noise with (^(t)) = and {^(t)^(t')) = 6{t — t'). Changing 
variables to s = — 1 + 2z, the coefficients for our case are 



a 



A{z) = ^ [-iV7s^ -2s + 7(2 + N) 



Biz) = - 



(l + s)(l-.) + -(l-7^) 



(40) 



The stochastic differential equation including noise is then obtained as 
1 ds ' 



a dt 



-N-fs-2s+-f{2+N) + ^ 



12 
a 



s) + -{l-is) 



e(t).(4i) 



A straightforward generalization to the multi-site case gives the following evolution 
equations 

1 dsi 
a dt 



-Ar7,(s)s2-2s, + 7^(s)(2 + Ar) 



+ 



12 

a 



[1 + Si){l - Si 



N 



1 - 7i(s)si 



(42) 



Neural networks using two- component Bose-Einstein condensates 



12 



The above equation is identical to (129 p up to the Gaussian noise term. This allows for 
the system to escape local minima in the energy landscape. The steady state evolutions 
then approach the correct thermal equilibrium averages as defined by the Boltzmann 
distribution. Combined with an annealing schedule, this may be used to find the ground 
state of the Ising Hamiltonian. As found previously, the due to the factors of N in (H2ll 
originating from bosonically enhanced cooling, the annealing rates may be made N 
times faster, allowing for an accelerated method of finding the solution to Ising model 
problems, as claimed in Ref. [5]. 

4. Learning and Pattern Completion 

In the context of the neural networks, learning and memory retrieval are typical tasks 
which illustrate their utility. Due to the equivalence of the BEC network to a continuous 
Hopfield model, this imphes that such processes should also apply to BEC networks. In 
this section we illustrate the equivalence between the two systems by simple examples 
of Hebbian learning and pattern completion. 

4-1. Learning 

The simplest example of unsupervised learning is the Hebbian learning rule |10]. Using 
the associations in Table [1] it is straightforward to write down the corresponding rule 
in the case of the BEC network. We follow the presentation give in Ref. [10] (sec. 2.5) 
for the case of continuous activation functions, since in the BEC system the measured 
spin is a continuous quantity. We assume that the BEC network starts with the system 
shown in Figure [1] with the Ising matrix set to 

Jij = 0. (43) 

The learning algorithm then proceeds as follows. We apply various magnetic field 

(n) 

configurations Xi where n labels the various pattern configurations that the network is 
exposed to during the learning process. Starting with the first field configuration n = 1, 
we apply this field and wait until the spins reach their equilibrium value, which will be 
given (for the first iteration) 

= <I>(/3AS^)) (44) 

We then update the Ising matrix according to 

J^j Jij + CSiXf' (45) 

where c is the learning constant which determies the speed of the learning process. We 
then make subsequent applications of the field \^^\ measure the field Si in each case, 
then make the replacement fHSj) . For n > 2, the spins do not simply take the values 
of dS]) since the fields Bi will in general be non-zero. The process is continued until 
the learning examples are exhausted, or the same set can be recycled. Other learning 
algorithms may be derived in a similar way using the associations in Table [1] 
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Figure 2. Pattern completion for a BEC network prepared with Hebbian learning for 
the two space invader characters in (a) and (b). The BEC network is evolved in time 
starting from the initial conditions as shown and evolve towards their steady state 
configurations. Parameters used are a — 1, N = 10^, fc^T = 0. 

4-2. Pattern completion 

As a simple example of the use of the BEC network for pattern completion, we test the 
set of equations given in (129|) using an Ising matrix trained using the Hebbian algorithm 
of the previous section. We numerically evolve a set of 13 x 10 equations forming a two 
dimensional grid at zero temperature from the initial conditions as shown in Figure |2j 
In Figures |2|l and |2b we start from fragments of the learned patterns while in Figure 
IJh we start from a randomly chosen spin configuration. We see that in all cases the 
spins evolve towards the learned configurations, with the BEC network completing the 
patterns as desired. For the random initial configuration, the spins evolve towards 
whichever configuration happens to be closer to the learned patterns. 

The time scaling behaviour is shown in Figure [3l In Figure |3t we plot the 
normalized Hamming distance 

^ = ^Ek^W-4"^l (46) 
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Figure 3. (a) Hamming distance D versus time for the pattern recognition task given 
in Figure[2^ for various boson numbers N as shown, (b) Times tg necessary for reaching 
a Hamming distance of e as a function of N . Dotted hne corresponds to 

between the evolved spin configuration and the learned spin configurations s^^\ We see 
that the general behavior is that the time for the pattern completion scales as ~ 
which can again be attributed to bosonic stimulated cooling. There is a logarithmic 
correction to this behavior, where there is an initial stiffness of the spins to move towards 
the steady state configuration. In Figure Wp we show the scaling of the time to reach a 
particular Hamming distance e with respect to N . We see that for large N all curves 
converge to the dominant ~ behaviour. This shows that BEC networks can equally 
well be used to perform tasks such as pattern completion, with the additional benefit of 
a reduced time in proportion to N . 

5. Summary and conclusions 

We have analyzed the BEC network proposed in Ref. [5] in terms of the theory of neural 
networks and found that it is equivalent to a stochastic continuous Hopfield model. In 
contrast to the continuous Hopfield model where the overall timescale of the evolution 
is determined by the capacitance within each unit, in the BEC network the timescale 
is determined via the rate of cooling. Due to bosonic stimulated cooling, the rate of 
cooling may be accelerated in proportion to the number of bosons on each site, which 
in turn accelerates the cooling rate of the entire system. The bosonic stimulated cooling 
makes the time evolution equations f l29|) on a single site not precisely the same as its 
Hopfield model counterpart ( fTSjl . but the difference merely gives a modification of the 
dynamics as the system heads towards equilibrium, the overall behaviour of the system 
as a whole remains the same. In particular, tasks such as pattern completion may be 
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performed using the BEC network, in the same way as the Hopfield model. 

In this context, it would appear that using a BEC network, rather than a physical 
implemetation of a Hopfield network, is nothing but a more complicated way of 
implementing what could be done equally well by either standard electronics or optical 
means |9]. Specifically, one could imagine using simply Hopfield circuits with small 
capacitances such that the timescale of the circuit is as small as desired. Other variations 
of optical implementations of Hopfield models allow for fast operation speeds. While for 
the zero temperature case this may be true, the BEC system does have the advantage 
that the random fluctuations following Boltzmann statistics is already built-in, and do 
not require additional circuitry to simulate. Another possible advantage is that the 
speedups can be made systematically faster by simply by increasing the number of 
bosons. A possible issue for the physical realization is whether the connections between 
each Ising site require response times of the order of the cooling time on each site. Apart 
from a simple slowdown due to bottlenecks in the transmission, such delay times in the 
information between each site can introduce instabilities in the system causing divergent 
behavior [T3]. We leave as future work whether the proposals in Refs. [7j and [8j can 
be treated with the same analysis. 
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